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ABSTRACT 


The  research  presented  in  this  report  is  both  a continuation  and  an 
extension  of  the  optimality  criteria  approach  to  structural  optimization 
reported  in  AFOSR-TR-7'3-1^51.  In  the  present  study  the  optimality  criteria 
method  is  extended  to  provide  a cai>ability  for  the  automated  minimum  weight 
design  of  elastic,  red-undEmt  structures  composed  of  one-  and  two-dimensional 
structural  elements  and  subjected  to  multiple,  Independent  static  loading 
conditions.  The  design  variables  are  teiken  to  be  the  thicknesses  of  the 
structural  elements.  These  variables  are  constrained  to  be  between  specified 
maximum  and  minimum  values,  as  are  the  internal  stresses  in  each  element  and 
the  nodal  displacements  of  the  stimcture.  Results  are  presented  to  Indicate 
both  the  excellent  performance  of  the  optimality  criteria  method  and  the  wide 
range  of  structures  which  can  be  designed  using  the  algorithm. 

Finally,  the  algorithm  is  extended  to  include  the  new  (to  automated 
design)  and  very  important  requirement  that  structural  integrity  under  the 
applied  loads  be  maintained  given  the  presence  of  existing  structursil  fatigue 
cracks.  This  requirement,  which  is  cast  in  an  energy  format  and  incorporated 
in  the  design  algorithm  as  an  inequality  constraint,  is  shown  to  have  a 
dominant  effect  in  the  design  of  safe,  mlnlmxim  weight  aircraft  structures. 
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The  research  presented  herein  is  berth  a continuation  and  an  extension 
of  the  optimality  criteria  method  developed  on  AI’'OSR  Grant  No. 
reported  in  AFOSR-TR- 75-1^51  and  subsequently  published  in  Ref.  1.  This 
research  activity  has  been  concerned  with  the  development  of  a capability 
for  automated  minimum  weight  design  of  elastic,  redundant  aircraft  struc- 
tures subjected  to  multiple  independent  static  loading  conditions  and  con- 
strained against  several  types  of  inadrlssable  structural  behavior.  Tne 
activity  reported  in  AFOSR-TR- 75-1^51  >^as  quite  fundamental  in  nature; 
the  basic  redesign  algorithm  was  developed  and  tested  on  elastic  truss- 
type  structures  subject  to  only  stress  constraints  and  constraints  on 
maximum  and  minimum  values  of  the  design  variables.  While  the  results 
were  very  encouraging  for  this  rather  restrictive  problem,  the  computational 
aspects  of  the  design  algorithm  were  not  completely  understood.  In  a 
broader  sense,  the  applicability  of  the  optimality  criteria  approach  to 
more  complex  (and  practical)  design  problems  was  untested.  Therefore,  in 
the  research  activity  which  followed  that  presented  in  Ai’OSR-TR-75-l^*^51> 
three  primary  objectives  have  been  to: 

(l)  Extend  and  refine  the  basic  oi>timality  criteria  method  to  improve 
its  already  efficient  performance  and  reduce  the  experience 
required  by  the  user  to  exercise  the  algorithm. 

(2  ) Increase  the  structural  modeling  capability  contained  in  the 
algorithm  to  the  point  where  two-  and  three-dimensional  assem- 
blages of  one-  and  two-dimensional  structural  elements  could  be 
designed  using  the  optimality  criteria  algorithm. 
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(3)  Extend  the  types  of  design  constraints  to  include  stress, 

displacement,  minimum  and  imiximura  constraints  on  the  design 
variables,  and  constraints  liased  on  gross  struct\U"aJ.  behavior, 
i.e.,  constraints  incluaed  in  other  existing  algorithms  for 
the  design  of  elastic  redundant  structures. 

To  a certain  extent,  these  three  research  objectives  were  interrelated 
since,  for  example,  algorithmic  performance  coxild  be  adequately  measured 
only  after  developing  design  experience  with  structures  more  ccmplex  than 
trusses  which  are  subjected  to  several  different  types  of  constraints. 

In  the  development  of  the  capability  for  treating  constraints 
against  inadmissible  structviral  behavior,  the  usual  constraints  on  stresses 
and  displacements  were  incorporated  into  the  design  algorithm,  ’lowever, 
since  a rather  substantial  literature  [2  - S]  demonstrates  the  ease  with 
which  optimality  criteria  methods  can  deal  with  constraints  on  gross  struc- 
tural behavior  (e.g.,  on  minimum  natural  frequencies  or  elastic  buckling 
loads)  a new  type  of  constraint  was  incorporated  into  the  algorithm;  a 
constraint  requiring  the  structure  to  maintain  its  integrity  given  the 
presence  of  existing  structural  fatigue  cracks.  This  fracture  mechanics 
constraint,  which  can  be  cast  in  a form  based  on  energy  concepts  for  the 
structure,  provides  a test  of  the  capability  of  the  algorithm  to  treat 
energy- type  constraints.  More  importantly,  the  development  of  this 
fracture  constraint  provides  a means  for  establishing  the  effect  that 
fracture  mechanics  considerations  can  play  in  determining  safe,  minimum 
weight  structural  designs.  As  will  be  shown  in  this  Report  , this  new 
fracture  mechanics  constraint  is  very  Important  and  in  fact  may  govern 


the  minimum  veight  structiu’al  design  or'  typical  aircraft  structural  sub- 
assemblies. 

This  Report  contains  a brief  review  of  the  mathematical  foundations 
of  the  optimality  criteria  method  utilized  in  this  research  in  Section  TI, 
as  well  as  a genereil  discussion  of  the  computational  features  incorporated 
in  the  design  algorithm.  Section  III  contains  a description  of  the  struc- 
tural modeling  capability  of  the  present  algorithm  and  also  a review  of  its 
computational  characteristics  for  several  well  known  test  problems.  The 
development  of  the  fracture  constraint  is  described  in  detail  in  Section  IV 
as  well  as  its  implementation  in  the  optimality  criteria  design  algorithm. 

A discussion  of  the  results  obtained  in  this  research  is  given  in  Section  V. 
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II.  THi;  NlATIiEMATICAL  OPTB-IALITY  CRITEHIA 


AND  REDESIGN  RULE 


A.  Revlev  of  the  Optltnallty  Criteria 


The  structural  weight  minimization  problem  is  stated  formally  as 

I 

follows:  find  a set  of  independent  design  variables  D , j = 1,2,...,N, 

J 

such  that 


^ = 1,2, ...,M 


(1) 


and 


W(D^) 


-5^  mm 


where  the  objective  function  V/(D.)  is  the  weight  and  the  constraint 

functions  g,  (D  ) are  behavioral  constraints  on  the  structural  response 
K j 

and  side  constraints  on  the  design  variables. 


The  necessary  criteria  for  optimum  design,  which  are  the  basis  for 
the  present  method,  are  the  Ktihn-Tucker  necessary  conditions  [6].  For  an 
optimum  Design  these  conditions  are. 


M 


\ IdT  = ^ = 1,2, ...,N 


(3) 


\ «k  = 


k = 1,2, ...,M 


(4) 


where  is  the  Lagrange  multiplier  associated  with  the  k-th  constraint. 
If  the  constraint  is  noncritical,  g^^  < 0 and  = 0;  if  the  constraint 
is  critical,  gj^  = 0 and  A^^  Thus,  the  Kuhn-Tucker  necessary 


conditions  may  be  simplified  to  read 


^ \ \ 

IF.  ^ / \ ~W.  ^ J 

J J 

where  ® 

The  active  side  constraints  on  the  design  variables  can  be  separated 
from  the  active  behaviorsil  constraints  and  the  Kuhn-Tucker  necessary  con- 
ditions for  a local  optimum  design  can  be  rewritten  as  follows : 


> 0 for  ft-11  K (<  M)  active  critical  constrainte. 


K 


Y 


k 


k 


0 


\ 


0 


\ 


0 


i ^ J, 


max 


i f- 


min 


(6) 

(7) 

(8) 


where  K is  the  number  of  active  behavioral  constraints,  J is  the  set 

of  design  variables  not  constrained  by  mlnijnum  or  maximum  side  values,  and 

J and  J . are  the  sets  of  design  variables  constrained  by  maxlnnim 
max  min 

and  minimum  values,  respectively. 


The  separated  form  of  the  Kuhn-Tucker  necessary  conditions  can  be 
rewritten  by  dividing  both  sides  of  Eqs.  (6)-  (8)  by  the  weight  gradient 
(assimied  positive,  the  case  for  most  structural  design  problems) 

tJ 

to  give 


K ^ 


B 

1, 

J 

€ 

J 

> 

1 , 

,1 

€ 

J 

max 

< 

1 , 

.1 

€ 

*^min 

(9) 


where  the  parameters  j - 1,2,...,N,  are  defined  as  the  design 


factors . 


The  relations  in  Eq.  (9)  are  the  criteria  for  minimiun  weight.  Tiiat 
is,  if  for  a given  design,  the  set  of  K active  behavioral  constraints  is 
such  that  the  correspondinf;  Lagrange  multipliers  are  positive  and  Eq.  >9) 
is  satisfied  for  sill  J,  then  the  design  satisfies  the  Kuhn-Tucker  neces- 
sary conditions  and  is  therefore  a local  optimum  design. 


P.  Redesign  Rule 


The  form  of  Kuhn-Tucker  necessary  conditions  given  by  Eq.  (9)  sug- 
gests an  associated  redesign  rtile.  Qualitatively,  the  criteria  of  Eq.  (9) 

imply  that:  (l)  if  T,  = 1.0,  then  D.  is  optimally  sized;  (?)  if 

J <3 

I,  > 1.0  and  D.  < D. 
j j j 


max 


then  D.  must  be  increased  to  obtain  an 

0 


improved  design;  and  (3)  if  I.  < 1.0  and  D.  > D.  , then  D must 

'^min  ^ 

be  decreased  to  obtain  an  improved  design.  Therefore,  the  iterative  re- 

(.X 

design  rule  is  taken  to  be  the  following:  given  the  o-th  design  1)  with 

K 

associated  design  factors  I.(D,  ),  then  the  design  variables  at  iteration 

J k 

(fj;  ^ 1)  are  related  to  the  design  variables  at  iteration  (.1  by  the 
relation 


J k J 

if 

•^min  ^ ^ 'max 

D. 

^max 

if 

^ '^max 

(10) 

•^min 

if 

''min 

where  ffl.o)  - 1,0.  Tn  the  present  study  the  function  f is  taken 
to  be 


r. 


* 


J k 


.1  K 


1/2 


= 

j 


where 


r 

J 


is  the  design  factor  at  ite>*Rtion  a. 


(11) 


This  redesign  rule  is  exact  for  statically  determinate  str'actures 
with  weight  functions  linearly  dependent  on  the  design  variables  and  sub- 
ject to  a single  frequency  or  buckling  constraint,  or  subject  to  a single 
load  condition  and  a single  stress  or  displacement  constraint.  The  rede- 
sign rule  is  approximate  for  the  design  of  statically  indeterminate  struc- 
tures and/or  multiple  load  conditions.  It  is  emphasized  that  the  design 
objective,  Eq.  (9),  is  exact. 


^ ! C . Conp'atational  Aspects 

1 

The  following  sequence  of  operations  defines  each  iteration  of  the 

I'  redesign  procedure:  (l)  identify  the  active  and  near-active  behavioral 

constraints  and  side  constraints;  (2)  calculate  the  gradients  and  l^a- 
gran^je  multipliers  of  the  active  behavioral  constraints;  (?)  evaluate  the 
design  factors  defined  by  Kq.  (9);  (^ ) either  terminate  the  redesign  pro- 

cedure If  Eq.  (9)  is  identically  satisfied  or  resize  the  structure  using  the 

; redesign  rule.  A brief  description  of  each  operation  is  described  below. 

■ 

Tlie  k-th  behavioral  constraint  is  strictly  active  if  g = O.U, 

K 

or  if  the  so-called  response  ratio,  K.  = g.  + 1.0  = 1.0.  However,  ^or 
numerical  purposes  a behavioral  constraint  is  identified  as  active  (or  near 
active)  and  governing  the  redesign  if  its  response  ratio  is  arbitrarily 
close  to  one  (l.o),  that  is,  if 


Li 


1.0-0  < < 1.0  (12  ) 

where  0 is  a preselected  constraint  buffer  parameter,  usually  contained 
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1 


in  the  range  O.OU  5^5  0*10*  identification  of  near-active  con- 

straints and  the  final  weight  obtained  are  sensitive  to  the  value  of  the 
constraint  buffer.  Large  values  of  p allow  early  Identification  of  near- 
active constraints  but  usually  result  in  heavy  final  designs  because  of  the 
large  spread  in  the  values  of  the  response  ratios  of  the  near-active  con- 
straints. Smaller  values  of  P usiislly  give  lighter  final  designs  becaiise 
of  the  comparitively  smaller  spread  in  near-active  constraint  response 
ratios . 

Analytic  expressions  for  the  components  of  the  constraint  gradients 
(that  is,  derivatives  of  stresses  and  displacements  with  respect  to  the 
design  variables)  are  obtained  by  differentiating  both  sides  of  the  stiff- 
ness relations  for  the  linear  elastic  structxire 


^K]ful  = fF] 


to  give 


^ ■ -I?  ■ '«)  • 


In  Eq.  (15),  [K]  is  the  system  stiffness  matrix,  [F]  is  the  vector  of 
equivalent  applied  nodal  forces  and  {u}  is  the  vector  of  developed 
nodeil  displacements.  The  apiplied  mechanical  loads  are  assumed  to  be 
independent  of  the  design  variables. 

Solving  Eq.  (l4)  for  the  displacement  gradient  gives 


^ - - IKl'^  ^ ful  . 


For  a given  design,  the  gradients  of  the  active  displacement  constraints, 

= (u^/u^)  - 1.0,  where  u^  is  the  i-th  component  of  the  displacement 


8 


vector  and  is  the  limiting  valxie,  is  obtained  by  reducing  the  product 

( >[K]/ ^Dj )ful  and  back  substituting  through  the  decomposed  stiffness 
matrix. 

The  gradient  of  active  stress  constraints  is  obtained  by  differenti- 
ating the  stress-displacement  relation  for  a member  to  give 

^ ^ 

J J 

Sind  substituting  for  the  displacement  gradient  from  Eq.  (l5).  In  P>q.  (l6), 
[a]  is  the  vector  of  element  stresses,  [D]  is  the  linear  elastic  con- 
stitutive matrix  and  [E]  is  the  linear  strain-displacement  matrix. 

The  calculation  of  the  Lagrange  multipliers  of  the  near-active  con- 
straints can  be  accomplished  in  several  ways.  In  the  present  study,  the 
Lagrange  multipliers  are  evaluated  by  minimizing  the  function  = ,?(a) 
where 


k .1 


The  conditions  for  minimizing  .5  are 


= 0,  k = 1,2, ...,K 


giving  a set  of  linear  algebraic  equations,  equal  in  number  to  the  number 
of  near-active  constraints,  to  solve  for  the  Lagrange  miiltipliers . 


Since,  in  general,  there  are  more  design  variables  than  active 
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constraints,  the  minimization  of  the  function  ,5  araovints  to  the  minimiza- 


tion of  the  mean  square  error  between  current  values  of  the  design  factors 


and  the  optlmTom  va].ue  l.C'  that  is, 


(l^  - l.O)'  . 


Move  limits  must  be  imposed  on  the  design  factors  before  resizing, 


because  the  redesign  rule  is  valid  only  near  the  optimum  and  is  approxi- 


mate elsewhere.  The  move  limits  are  taken  as  uniform  so  that 


1.0-  |a|  < I < 1.0+  |Ai,  J=1,2,...,N 

J 


where  |a|  is  a preselected  move  limit  parameter. 


Resizing  the  structure  may  cause  the  set  of  current  active  re- 


straints to  change  (constraint  switching);  excessive  constraint  switching 


may  cause  infeasible  designs.  Omall  move  limits  minimize  this  effect. 


However,  constraint  switching  and  infeasible  designs  may  still  occur. 


particularly  in  the  initial  stages  of  the  redesign  where  the  resizing 


rule  is  most  approximate.  Feasible  designs  are  then  obtained  by  scaling 


the  design  variables  so  that 


. E 


where  R is  the  response  ratio  of  the  most  active  (and  violated) 


constraint. 


For  problems  with  stress  limits  only  the  move  limits  must  be  rela- 


tively small  with  O.OU  5 1^1  0.08.  (in  addition,  for  problems  with 


stress  limits  only,  a fully  stressing  move  is  made  on  those  design  varia- 


bles with  corresponding  design  factors  less  than  zero.  This  move  is 


useful  in  quickly  Identifying  the  set  of  design  variables  constrained  by 
minimum  values.)  Kor  displacement  - limited  problems  the  move  limits  may 
be  greater  and  are  typically  |a|  =0.25.  fhe  set  of  near-active  con- 
straints in  displacement-limited  problems  is  usually  smaller  and  more 
stable  than  the  set  in  stress-limited  problems;  therefore,  relatively 


t'V 

R. 


larger  move  limits  are  possible  in  displacement-limited  design. 


Computational  experience  indicates  that  the  following  default 
values  for  the  two  preselected  design  parameters  P and  |a|  are 
usually  adequate  for  obtaining  near  opt.lmum  designs;  (l)  for  stress- 
limited  problems,  p = O.OB,  |a|  = 0.08;  (2)  for  displacement-limited 

problems  (and  stress  end  displacement  limited  problems),  (i  = 0.10,  |5]  = 

0.25.  It  is  usually  possible  to  improve  the  final,  designs  obtained  with 
the  defaiilt  parameters  by  reducing  the  value  of  the  constraint  buffer  P, 
although  more  iterations  may  be  necessary  for  convergence.  For  displacement- 
limited  problems  the  buffer  P can  be  chosen  independently  of  the  move 
limit  |a1.  However,  because  of  the  sensitivity  of  the  stress  constraints 
in  stress-limited  problems,  it  is  necessary  to  reduce  the  move  limit 
parameter  as  well  as  the  constraint  buffer  parameter. 


Finally,  it  is  necessary  to  modify  the  functional  .9  and  recalculate 
the  Lagrange  multipliers  and  the  design  factors  if,  in  the  trial  resizing, 
(1)  new  side  constraints  are  identified,  or  (2)  negative  multipliers  are 
obtained.  When  a new  side  constraint  is  identified,  it  is  necessary  to 
delete  the  design  variable  from  the  sum  in  Eq.  (l7)  over  the  variables 
contained  in  the  set  J.  When  a negative  multiplier  is  obtained,  it  is 
necessary  to  delete  the  corresponding  constraint  from  the  sum  in  Eq.  (17) 
over  the  set  of  K active  behavioral  constraints. 
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Both  the  deletion  of  design  variables  associated  vd.th  side  constraints 
and  the  deletion  of  behavioral  constraints  associated  with  negative  multi- 
pliers may  require,  within  each  redesign  cycle,  one  or  more  iterations  (one 
complete  solution  per  iteration)  to  obtain  a stable  list  of  side  constraints 
and  active  behavioral  constraints.  When  the  constraint  lists  stabilize, 
the  resized  structure  is  analyzed  and  scaled  and  the  redesign  procedure 
reentered . 

D.  Algorithm  Efficiency 

The  optimality  criteria  method  of  the  present  study  is  in  general 
less  efficient  than  previous  physical  optimality  criteria  methods.  How- 
ever, the  present  method  is  more  general  than  previous  optimality  criteria 
methods  since;  (l)  convergence  to  local  optimum  designs  is  obtained;  and 
{?  ) special  forms  for  different  behavioral  constraint  types  eire  not  needed. 
The  generality  of  the  present  method  is  obtained  at  a greater  computational 
cost,  primarily  for  constraint  gradient  and  Lagrange  mviltiplier 
calculations. 

The  present  method  is  more  efficient  than  mathematical  programming 
methods,  primarily  because  of  the  relative  simplicity  of  the  resizing 
calciilation.  The  efficiency  advantage  is  most  prevalent  for  design 
problems  in  v;hich  the  number  of  final  active  constraints  is  less  than  the 
niaraber  of  design  variables.  However,  the  results  for  the  65-bar  truss 
stress-limited  problem  [7l  show  the  present  method  with  a significant 
(2  to  l)  efficiency  advantage  over  the  latest  generation  mathematical 
programming  methods  for  problems  with  a large  number  of  design  variables 
and  a large  number  of  final  active  constraints. 
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The  computer  code  for  the  present  method  incorporates  a number  of 


•! 

1 


features  to  reduce  both  computer  storage  requirements  and  execution  time. 
Briefly,  thesp  features  are ; 

(l)  dynamic  allocation  of  all  arrays  in  terms  of  the  relevant 
control  parameters; 

(B  ) crvervriting  of  arrays  ; 

(3)  storage  of  only  non-zero  degrees  of  freedom; 

(it)  design  variable  linking; 

(3)  system  stiffness  calculation  by  stiffness  gradients  with 


[K]  = 


where  the  stiffness  gradients  are  calc\ilated  only  once  and 
stored  in  compacted  form; 

(6)  analytic  gradient  calculations;  and 

(7)  elimination  of  dependent  gradients  (due  to  symmetry  related 
design  variable  linking). 

Additional  details  are  given  in  the  Appendix  A. 


Further  and  significant  efficiency  gains  could  be  obtained  by: 

(1)  incorporating  approximation  concepts  to  reduce  the  number  of  complete 
structural  analyses  (tiiat  is,  decomposition  of  the  system  stiffness  matrix 
and  reduction  and  back  substitution  of  the  load  vectors);  and  (7)  using 


iterative  solution  methofis  to  reduce  the  ntunber  of  complete  solutions  for 


the  Lagrange  miHtipliers . 
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E.  Comparison  with  Other  Optimality  Criteria  Method; 


Optimality  criteria  methods  using  the  Kuhn-Tucker  necessary  condi- 
tions as  the  criteria  for  optimum  design  have  been  proposed  [4].  However, 
the  associated  redesign  rules  in  these  methods  do  not,  in  general,  con- 
verge to  designs  satisfying  the  Kuhn-Tucker  necessary  conditions  because: 
(l)  elements  governed  by  stress  constraints  are  resized  usir^  stress 
ratio;  and  (2  ) in  most  cases  the  "envelope  method"  is  used  wherein  re- 
sizing is  accomplished  considering  the  most  critical  constraint  for  each 
element . 

Optimality  criteria  methods  with  redesign  rules  capable  of  converg- 
ing to  designs  satisfying  the  Kuhn-Tucker  necessary  conditions  were  pro- 
posed by  Kiusalaas  [3]  (and  implemented  for  multiple  constraints  by  Rizzi 
[8]  ) and  by  Teral  [9l«  The  method  of  Kiusalaas  differs  from  the  method 
of  the  present  study  primarily  in  the  calculation  of  the  Lagrange  multi- 
pliers of  the  active  constraints.  In  Ref.  8 a set  of  simultaneous  linear 
equations  for  the  Lagrange  multipliers  is  obtained  by  requiring  strict 
activity  of  constraints  identified  as  near-active  during  resizing;  that 
is,  g^  = ~ 0.0.  The  set  of  equations  is  solved  using 

f)auss-Seidel  iteration.  Inactive  constraints  corresponding  to  negative 
Lagrange  multipliers  are  deleted  when  encountered  and  the  side  constraint 
list  is  updated  before  resizing. 

The  results  for  design  problems  with  a small  and  stable  list  of 
final  constraints  shows  the  method  of  Ref.  8 to  converge  to  designs  satis- 
fying the  Kuhn-Tucker  necessary  conditions  with  all  final  constraints 
strictly  active  (and,  therefore,  capable  of  giving  final  designs  lighter 
in  weight  than  final  designs  obtained  using  the  present  method).  However, 


the  method  of  Ref.  8 has  no  advantage  over  the  method  of  the  present  study 
for  problems  having  a set  of  active  constraints  which  is  neither  small  nor 
stable  as,  for  example,  in  large  primarily  stress-limited  problems. 

Indeed,  the  potential  efficiency  advantage  of  optimality  criteria  methods 
over  mathematical  programming  methods  may  be  lost  when  applying  the  method 
of  Ref.  8 to  large  problems.  The  results  of  Ref.  8 imply  that  infeasible 
designs  with  a large  number  of  active  and  violated  constraints  are  needed 
from  the  beginning  to  properly  sort  the  final  constraint  set.  This  re- 
quires the  solution  of  large  sets  of  linear  equations  for  the  determination 
of  the  Lagrange  multipliers  throughout  the  design  process. 
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III.  EXTENDED  MODELING  CAPABILITY 

The  generality  of  the  optimality  criteria  method  of  the  present  study 
is  demonstrated  by  resiilts  for  the  design  of  trusses,  plates  and  stiffened 
plate  systems.  These  results  were  obtained  using  a computer  program  coded 
in  FORTRAN  IV  for  use  on  the  UCI.A  L'ampus  Computing  Network  IBM  36o/91 
system.  All  calculations  were  performed  in  a single  precision  arithmetic; 
the  FORTRAN  G compiler  option  was  used,  unless  noted  otherwise.  The  progreim 
structure  is  described  in  Appendix  A and  required  data  cards  are  listed  in 
Appendix  B. 

A . Finite  Element  Library 

Five  finite  element  types  are  included  in  the  present  computer  program 
for  n\imerical  stress  analysis.  These  element  types  are: 

1.  the  two-node  truss  element; 

2,  the  three-node  subparametric  truss  element; 

5.  the  planar  constant  strain  triangle  (CST); 

A.  the  planar  four-node  isoparametric  quadrilateral  (QUADI); 

and 

5.  the  planar  eight-node  isoparametric  quadrilateral  (QUAD8). 

i 

The  more  complex  isoparametric  elements  are  preferred  in  nearly  all  I 

’ 

cases  over  the  simple  constant  strain  triangle  for  planar  analysis.  The 
advantages  of  the  isoparametric  elements  are : 

1.  improved  structural  modeling  capability; 

2..  improved  accuracy  over  simple  triangles  and  rectangles  for  the 


The  primary  disadvantage  of  the  complex  elements  is  the  extra  effort 
required  to  generate  the  element  stiffness  matrix.  However,  since  element 
stiffness  matrices  are  calculated  only  once  (to  form  the  stiffness  grad- 
ients), the  additional  cost  of  the  initial  stiffness  calculation  is  not 
significant  in  a design  program. 


Element  stiffness  matrices  are  obtained  by  mimerical  integration. 

The  complex  elements  are  mapped  into  simple  shapes  in  a natural  coordinate 
system  in  which  the  stiffness  coefficients  are  evaluated  using  Gaussian 
quadrature.  (Two-point  quadrature  is  used  for  the  three- node  truss  ele- 
ment, two-  and  three-point  formulas  are  used  for  the  QUAlA,  and  two-, 
three-  and  four- point  formulas  are  used  for  the  QUAD8  element.)  Additional 
details  on  stiffness  generation  are  included  in  Appendix  C. 

Note  that  it  is  not  necessary  to  retain  a constant  strain  triangle 
(as  is  done  in  the  present  program).  A straight-sided  triangle  may  be 
formed  by  collapsing  one  side  and  two  nodes  of  the  QUA.oU  to  a single  node; 
a curved-sided  triangle  may  be  formed  by  collapsing  one  side  and  three  nodes 
of  a QUAI)8  element  to  a single  node.  The  CST  element  is  presently  used 
only  in  the  transition  zone  between  fine  and  coarse  element  grids. 


Four  behavioral  constraint  types  are  included  in  the  present  computer 
program.  These  constraint  types  are: 

1.  limits  on  allowable  element  stresses; 

?■ . limits  on  allowable  nodal  displacements; 

5.  limits  on  allowable  relative  nodal  displacements;  and 
J(.  fracture  mechanics  constraints,  i.e.,  limits  on  developed  open- 
ing mode  stress  intensity  factors. 
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Application  of  mathematical  optimality  criteria  methods  to  bucKlin*^ 
and  frequency  constraints  has  l^een  demonstrated  i)y  Ternl  I'-)]  (with  special 


forms  of  the  redesign  rule  for  problems  with  oaly  one  active  constraint 
type).  Application  to  aeroelastic  flutter  constraints  has  been  demon- 
strated  by  Rizzi  [8]  so  that  these  results  are  not  duplicated  in  the 
present  study.  Rather,  a new  constraint  type,  against  fast  fracture,  is 
considered.  The  need  to  include  fractiire  constraints  in  analysis  aiid  de- 

2 f 

sign  is  evidenced  by  the  increasing  volume  of  literature  on  the  subject  in  ; 

' - t 

a number  of  structural  engineering  problems,  ranging  from  nuclear  reactor  1/ 

' components  to  lightweight  aircraft  components.  Little  work  has  been  done 

to  rationally  include  fracture  constraints  in  automated  design.  This  import- 
ant design  consideration  is  described  in  detail  in  Section  IV. 

B.  Truss  Results 

Both  two-  and  three-dimensional  truss  problems  were  investigated  using 
the  present  methods.  Results  for  the  10,  63  and  7Z  bar  trusses  with 

stress  limits  only  are  given  in  Ref . 7.  Results  for  the  10,  and  63  bar 
trusses  with  stress  and  displacement  limits  are  given  in  Ref . 1.  These  re- 
sults compare  favorable  with  previous  results  obtained  by  both  mathematical 
programming  and  physical  optimality  criteria  methods.  The  present  method 
is  effective  in  generating  near  optimum  designs  and  is  capable  of  converg- 
ing to  true  optimum  designs  when  the  physical  optimality  criteria  methods 
cannot.  The  present  method  is  also  very  efficient  compared  to  the  most 
recent  design  methods  based  on  mathematical  programming.  Hie  number  of 
complete  structural  analyses  needed  in  the  present  method  are  essentially 
the  sam»^  as  the  numl-er  needed  in  the  mathematical  programming  methods,  but 
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the  redesign  procedure  is  much  less  costly.  Tiiis  efficiency  is  retained 
for  large  problems  as  illustrated  by  the  results  for  the  5-bar  truss 
wing  carry  through  box.  Tliese  results  are  reviewed  because  of  the  Imp^re- 
ance  of  this  problem  in  the  literature. 


6j5-bar  tniss  wing  carry  through  box,  two  load  condltiona, 
stress  limits  only. 

Tlie  structural  configuration  is  shown  in  Fig.  1.  Hie  material  is 
titanium  with  E = l6.0  x 10^  psi.,  o = 0.l6  pci.^  and  = -r  100. 0 x IC^ 
psi.  No  design  variable  J.inking  is  specified;  therefore,  this  pro>lem  has 

p 

G’j  independent  design  variables.  A minimum  area  constraint  - 0.01  in.  is 
imposed  on  the  design  variables;  maximum  limits  are  not  imposed. 

The  problem  was  initiated  from  a uniform  starting  design.  A weight 
ol'  51^*5  Ihs.  was  obtained  after  15  analyses  using  the  stress-limited  n 

default  design  i^arameters  P = O.OB  and  |a1  = 0.08.  An  improved  res'ult 
with  W = 50'’^^  lbs.  was  obtained  after  lx-  analyses  with  & = O.OP  and 
|a|  = 0.05.  These  results,  detailed  in  Table  1,  car.pare  favorably  with  the 
stress  ratio  design  of  505^+  lbs.  after  50  analyses  and  the  design 
obtained  usirig  mathematical  programming  llO]  with  W = 4o?7  lbs.  after 
14  analyses. 

It  is  clear  from  Table  ?.  that  the  convergence  of  the  present  method 

is  not  monotonic.  The  fluctuation  in  the  number  of  active  constraints  is 

primarily  due  to:  (l)  the  approximations  inherent  in  the  redesign  rule; 

{?.  ) the  sensitivity  of  member  stresses  to  changes  in  member  areas;  and 
(3)  the  large  number  of  active  constraints  in  the  final  design.  Hie 
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final  designs  in  Refs.  ^ and  10  are  essentially  fully  stressed.  The  final 
design  in  Ref.  10  has  active  stress  constraints  and  five  active  minimurr. 

area  constraints.  The  final  design  of  the  present  study  has  active 

stress  constraints  and  nine  active  minimum  area  constraints;  this  design 
is  not  completely  converged  since  the  Kuhn-Tucker  necessary  conditions  are 
not  identically  satisfied.  However,  the  square  root  of  the  design  factors 
of  UR  of  the  5^^  nonminim\im  area  bars  are  converged  to  + 5 percent  of 
unity  and  the  remaining  design  factors  are  converged  to  + 10  percent  of 
unity.  Continued  redesign  does  not  lead  to  large  changes  in  variable  sizes 
or  final  weight. 

The  final  results  for  this  problem  using  the  present  method  required 
l8  CRJ  sec.  for  execution  on  the  IBM  560/91,  using  the  FORTRAIJ  H compiler. 
The  method  of  Ref.  10  required  approximately  48 CPU  sec.  of  execution  time 
on  the  same  IBM  560/91  computer  and  the  FORTRAN  H compiler  for  a one  percent 
overweight  design  (V/  = 5045  lbs.),  using  default  design  parameters.  Subse- 
quent results  in  Ref.  10  with  design  parameters  selected  to  speed  converg- 
ence required  approximately  56  CRJ  sec.  for  a one  percent  overweight 
design  (W  = 50R5  lbs.).  Thus,  the  present  optimality  criteria  method 
retains  a computational  efficiency  advantage  over  m.athemat ical  programming 
methods  even  when  the  number  of  design  variables  is  larger  than  the  number 
of  degrees  of  freedom  and  the  niimber  of  final  active  constraints  is  not 
small  compared  to  the  number  of  design  variables. 

C.  Center  Hole  Plate  Results 

The  rainljirum  weight  design  of  a center  hole  rectangular  plate  subject 
to  a uniform  tensile  stress  is  considered.  The  structural  configuration 
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is  shovm  In  Fig.  ?. . Three  design  cases  are  considered:  (l)  xinreinforced 

hole  with  stress  and  minimum  side  limits;  (;: ) unreinforced  hole  with 
stress,  minimum  and  maximum  side  limits;  and  (5)  reinforced  hole  with 
stress  and  minimum  side  limits.  RLane  stress  conditions  are  assumed. 

Due  to  symmetry,  only  a quarter  of  the  plate  need  be  considered.  The 
finite  element  Idealization  of  the  quarter  plate  is  shown  in  Fig.  3.  The 
skin  is  modeled  by  29  planar  four-node  isoparametric  quadrilaterals  and 
the  reinforcement  by  five  two-node  trtiss  elements.  No  design  variable 
linking  is  imposed  and  there  are  5^  design  variables,  the  thicknesses  of 
the  29  planar  elements  and  the  areas  of  the  five  truss  elements.  Numeri- 
cal integration  for  the  stiffness  is  accomplished  by  two-point  Gaussian 
quadrature.  The  largest  of  the  integration  point  Von  Mises  equivalent 
stresses  of  an  element  is  assumed  to  govern  the  design  of  the  element. 
Initial  stress  analysis  of  the  unreinforced  plate  is  checked  against  re- 
sults obtained  using  the  finite  element  model  with  I8  eight-node  quadri- 
laterial  elements  (Fig.  k). 

For  an  three  cases,  the  applied  tensile  stress  is  cr  = 21,70U.t 
psi.  and  the  material  is  alviminum  7075  (with  E = 10.8  x 10^^  psi.,v  = 
0.335,  n = 0.1  pci.  and  ~ ^ psi.).  A minimum  side  limit  = 

0.01  in.  is  iinposed  in  all  cases.  All  results  were  obtained  using  the 
stress  limited  defa\ilt  parameters  P = 0.08  and  |a|  = O.08.  All  results 
appear  to  compare  favorably  with  the  design  presented  by  Pope  (U],  but 
only  a qualitative  comparison  is  possible  since  Pope  does  not  provide  load- 
ing or  material  property  information.  Pope's  study  also  was  slightly  dif- 
ferent since  his  design  variables  were  the  nodal  thicknesses. 

2I+ 
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C.l  IJnreinforced  hole,  stress  and  mlnlnnm  limits 
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The  problem  was  initiated  from  a uniform  design  with  skin  thickness  = 
0.084  in.  A weight  of  0.5l8  lbs.  was  obtained  in  l8  iterations.  The 
resiilts  given  in  Table  5 shcrw  the  large  buildup  in  thickness  around  the 
hole. 

The  design  is  not  fully  converged.  There  are  20  active  stress  con- 
straints and  five  active  minimum  side  constraints.  The  square  root  of  the 
design  factors  ' are  converged  to  1.0  ♦ 5 X>ercent  (with  the  largest 
variation  in  the  design  variables  close  to,  but  not  equal  to,  the  mlninimi 
area).  Therefore,  as  in  the  previous  examples,  continued  redesign  would 
not  lead  to  large  changes  in  the  design  variables  or  final  weight. 

C.2  Unreinforced  hole,  stress,  minimum  sind  maximum  limits 

The  large  thickness  buildup  around  the  unreinforced  holes  is  unreal- 
istic for  plate  manufacture  and  may  also  violate  the  conditions  for  plane 
stress  analysis.  These  difficulties  can  be  minimized  by  imposing  raaxim’jm 
as  well  as  minimum  limits  on  the  design  variables.  The  procedure  was 
initiated  from  a uniform  starting  design  with  thickness  = 0.084  in.  The 
design  problem  was  modified  by  imposing  a maximum  side  limit  = 0.12  in. 
on  the  design  variables.  A weight  of  0.325  lbs.  was  obtained  in  19 
iterations  with  final  results  as  given  in  Table  4.  This  weight  is  as  ex- 
pected greater  than  the  design  obtained  with  no  maximum  side  limits. 

The  design  is  fully  converged  with  21  active  stress  constraints, 
one  maximum  side  constraint  and  six  minimum  side  constraints.  Hie  design 
factors  for  the  elements  not  constrained  by  minimum  or  maximum  side  limits 
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MINIMUM  WEIGHT  DESIGN  OF  PLATE  WITH  UNREINFORCED 
CENTER  HOLE,  STRESS  AND  MINIMUM 
THICKNESS  LIMITS 


! Variable 
1 

Thickness 

Variable 

Tiiickness 

T ^ 

1 

C.1R7O 

18 

0.0162 

i 

0.106't 

17 

0.0166  ^ 

5 

0.05‘^R 

18 

0.0230  ' 

L 

0 .0100 

19 

0.0464 

•3 

0.0100 

20 

0.0100 

6 

0.01 1+1 

21 

0.0141 

7 

0.071+9 

2? 

0.0182 

R 

0.07R9 

25 

0.0210 

9 

0.0524 

24 

0.0235 

10 

0.0100 

25 

0.0550  j 

11 

0.0107 

26 

0.0186  1 

IL’ 

0.0175 

27 

0.0251  j 

0.0582 

28 

0.0286  ' 

1 

Ih 

0.0459 

29 

0.0207  ' 

1 

^ 1 

0 .0100 



v;  = 0.3179  lbs.  j 
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TABU'!  4 

MINIMLFM  TOIGHT  DESIGN  OF  PLATE  WITH  UNRi;iNFORCED 

CENTI':K  hole,  stress,  minimum  and  maximum 
T>{ICKNESS  LIMITS 


* 


r-N 

1.  w 


I 


I 


are  identically  equal  to  1.0,  the  design  factors  for  the  element  con- 
strained hy  the  maximinn  limits  is  greater  than  ''  .0  and  the  design  factor 
for  the  elements  constrained  by  minimum  side  limits  are  less  than  1.0. 

C . 3 Reinforced  hole,  stress  and  minimum  limits 

The  unrealistic  thickness  buildup  near  the  hole  of  the  iinreinforced 

plate  may  also  be  averted  by  reinforcing  the  hole.  The  reinforcement  is 

modeled  by  five  two-node  truss  elements.  The  design  process  was  initiated 

from  a starting  design  with  uniform  thickness  = in.  and  uniform  re- 

2 

inforcement  area  =0.36  in  . A weight  of  0.3l6  lbs.  was  obtained  in  13 
iterations;  details  of  this  design  are  given  in  Table  3.  As  e>q^cted,  the 
presence  of  the  reinforcement  causes  a reduction  in  the  skin  shlckness 
around  the  hole  and  the  final  weight>  as  compared  to  +he  unreinforced  case. 
(However,  in  all  cases,  the  difference  in  the  design  variables  and  final 
weights  is  not  large. ) 

This  design  is  not  fully  converged.  For  this  problem  there  are  ?.h 
active  stress  constraints  and  seven  active  minimum  constraints.  The 
square  root  of  the  design  factors  are  converged  to  1.0  + 2 percent. 
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TABLE  5 


MINIMUM  WEIGHT  DESIGN  OF  PLATE  WITH  REINFORCED  CENTER  HOLE 
STRESS  AND  MINIMUM  THICKNESS  LIMITS 


Variable 

Thickneafi 

Variable 

niickness 

1 

0.0802 

.19 

0.0387 

O 

C. 

O.L225 

20 

0.0100 

\ 5 

0.0i*l8 

21 

0.0158 

I 

0.012.2 

22 

0.0178 

5 

0.0100 

23 

0.0'^  17 

6 

0.0241 

24 

0.0246 

7 

0.0453 

25 

0.0262 

8 

0.0867 

26 

o.cr2o6 

9 

0.0356 

27 

0.0258 

10 

0.0100 

28 

0.0262 

11 

0.0168 

29 

0.0229 

■ 

12 

0.0232 

Area 

13 

0.0352 

30 

0.0717 

li+ 

0.0618 

31 

0.0182 

1*3 

■■^.0100 

32 

C.OlOO 

16 

0.018s 

33 

'.\0100 

17 

0.0195 

34 

0.0100 

I 


IV.  OPTIMIZATION  INCLUDING  FFACTliPJC  CONSTRAINTS 


The  need  to  include  fracture  constraints  in  structural  analysis  end 
design  is  evidenced  by  the  increasing  volime  of  literature  on  the  subject 
in  a number  of  structural  engineering  problems,  ranging  from  nuclear  re- 
actor components  fli’ - 15]  to  lightweight  aircraft  components  [l6-?0]. 
iiowever,  little  work  has  been  done  to  rationally  include  fracture  con- 
straints in  automated  structural  design.  Simultaneous  failure  concepts 
were  used  by  Hxint  [l^]  to  design  stiffened  panels.  H'lnt  assumed  an  opti- 
mum minimum  weight  design  would  occur  when  the  stress  in  the  frame /crack 
stopper  and  the  average  stress  in  the  skin  panel  containing  the  crack  sLm- 
liltaneously  reached  allowable  values.  Recently,  !)avis  ' ll  used  matnema- 
tical  programming  methods  to  design  an  integrally  stiffened  two-panel  tox 
beam  with  stress,  displacement,  buckling  and  fracture  constraints,  as  well 
as  side  constraints  on  the  design  variables.  The  design  variables  were 
skin  thickness  and  the  thickness,  depth  and  spacing  of  the  integral  stiff- 
eners. The  fracture  constraints  included  limits  on  developed  stress  intens- 
ity factors  under  large  static  loads  as  well  as  limits  on  crack  grovrt.h  dur- 
ing cyclic  loads.  Analytical  results  were  used  for  the  stress  intensity 
factor  under  static  loads  and  for  the  increment  in  the  stress  intensity 
factor  under  cyclic  loads. 

The  simultaneous  failure  method  of  Ref.  l6  has  limited  application 
since  simultaneous  failure  design  may  not  yield  a true  optimum  design  as, 
for  example,  when  side  constraints  on  the  design  variables  are  included. 

The  mathematical  programming  method  of  Ref.  R1  is  capable  of  converging  to 
true  optimum  designs,  but  the  method  is  computationally  expensive,  even 


when  implemented  with  analytical  resxilts  for  developed  stress  intensity 
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factors.  Extension  of  the  method  of  Ref.  Ml  to  problems  requiring  n\imeri- 
oal  calculation  of  developed  stress  intensity  factors  would  be  computation- 


ally prohibitive. 

Tlie  pi-esent  optimality  criteria  method  can  be  used  for  efficient 
automated  minimum  weight  design  of  structural  components  for  which 
analytic  solutions  for  developed  stress  intensity  factors  are  not  avail- 
able. The  finite  element  method  is  used  for  stress  analysis  and  the  strain 
energy  release  rate  method  (the  compliance  method)  is  used  to  calculate 
developed  opening  mode  stress  intensity  factors.  Only  two  structixral 
analyses  are  needed  at  each  design  iteration  to  calculate  necessary  re- 
sponse gradient  information  as  well  as  the  developed  stress  intensity 
factor.  The  method  is  applicable  to  design  \inder  large  static  loads,  includ- 
inf  fail-safe  design  of  aircraft  fuselage  panels. 

A brief  review  of  the  principles  of  linear  elastic  fracture  mechan- 
ics (LEFM)  follows,  after  vrhich  is  presented  the  response  gradient  calcu- 
lations required  to  incorporate  the  fracture  mechanics  constraints  into  the 
design  algorithm.  Numerical  results  for  a fuselage  design  problem  are 
used  to  illustrate  the  importance  of  fracture  mechanics  constraints. 

A . Irwin  Theory  of  Linear  tHastlc  Fracture  Mechanics 

Irwin  [2?]  has  shown  that  a general  linear  elastic  stress  state 
asymptotically  close  to  a crack  tip  is  a linear  combination  of  three 
distinct  stress  states,  each  stress  state  associated  with  a partic’ilar  mode 
of  crack  tip  deformation.  The  three  deformation  modes  are  Mode  I,  the 
opening  mode  in  which  the  crack  surfaces  separate  in  a direction  normal  to 
the  crack  surface;  Mode  II,  the  sliding  mode  in  which  the  crack  surfaces 
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slide  along  a direction  perpendicular  to  the  crack  front;  and,  M.ode  IT,  tne 
tearing  mode  in  which  the  crack  surfaces  slide  in  a direction  parallel  to 
the  crack  front. 


The  near  crack  tip  stresses  and  displacements  in  each  mode  have  been 
obtained  using  complex  variable  techniques  [23]  and  eigenfunction  expansions 
[2U].  The  near  crack  tip  results  are  the  product  of  an  intensity  factor 
and  functions  of  the  position  coordinates.  The  results  for  Mode  I deforma- 
tion of  an  infinitely  sharp  crack  are 
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where  is  the  opening  mo<ie  stress  intensity  factor  and  x and  y 

are  cartesian  coordinates  and  r and  S are  polar  coordinates  v/ith 
origin  at  the  crack  tip  (Fig.  3).  In  Eq.  (2U),  u and  v are  the  dis- 
placements along  X and  y,  respectively,  u is  the  elastic  shear 
modulus  and 
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plane  strain 

generalized  plane  stress 
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Tliese  restilts  are  valid  only  for  r « c where  c is  the  crack  half- 
length  . 

The  stress  intensity  factor  controls  the  intensity  of  the  elas- 

tic stresses  near  the  crack  tip.  For  a small  extent  of  yielding  at  the 
crack  tip,  the  elastic  stresses  outside  the  plastic  zone  are  not  greatly 
disturbed  and  are  completely  described  by  the  stress  Intensity  factor. 
Therefore,  the  stress  intensity  factor  is  the  relevant  fracture  parameter 
for  elastic  and  limited  yielding  fracture. 

The  value  of  the  stress  intensity  factor  is  dependent  on  the  load- 
ing and  specimen  geometry  and  is  obtained  by  an  analjd;ical  or  numerical 
solution  of  the  boundary  value  problem.  Thus,  the  near  crack  tip  stresses 
are  the  same  for  different  boundary  conditions  if  the  stress  Intensity 
factors  are  the  same. 

The  criterion  for  the  onset  of  fracture  is  that  the  state  of  stress 
at  the  crack  tip  reaches  a critical  state,  or  equivalently,  that  the  stress 
intensity  factor  reaches  a critical  value.  Therefore,  for  the  opening  mode, 
fast  fracture  first  occurs  when 


(26) 


where  is  the  critical  opening  mode  stress  intensity  factor. 

Experimental  results  show  the  critical  opening  mode  stress  intensity 
factor  to  be  primarily  dependent  on  specimen  thickness  and  the  state  of 
stress  at  the  crack  tip.  For  thick  test  specimens,  predominantly  in  a 
state  of  plane  strain,  the  critical  opening  mode  stress  intensity  factor 
is  a material  constant  and  is  denoted  The  critical  plane  strain 
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stress  Intensity  factor  can  be  measured  for  one  carefully  controlled  loadin*i 
and  geometry.  The  corresponding  critical  state  of  near-crack-tip  stresses 
is  the  same  for  all  other  cases  of  loading  and  geometry. 

For  test  specimens  not  predominantly  in  a state  of  plane  strain  the 
critical  opening  mode  stress  intensity  factor  is  not  constant,  but  increases 
with  decreasing  thickness.  The  critical  state  of  stress  is  not  independent 
of  the  boundary  conditions.  Therefor,,  the  test  specimen  and  the  service 
specimen  must  be  of  similar  geometry  and  size. 

R.  Solution  for  Developed  Stress  Intensity  Factors 

A number  of  methods  can  be  used  to  solve  for  developed  stress  intens- 
ity factors  (given  the  finite  element  solution  of  the  boundary  value  prob- 
lem. These  methods  include : 

(l)  stress  matching; 

(2  ) displacement  matching ; 

(3)  energy  release;  and 

(^)  embedded  singularity. 

These  methods  are  briefly  reviewed  with  an  emphasis  on  the  calcula- 
tion of  the  developed  opening  mode  stress  Intensity  factor  K^. 

B.l  Stress  matching 

In  this  method,  corresponding  values  of  the  stresses  and  position 
coordinates  are  substituted  into  the  singular  solution,  Eq.  (23),  to  solve 
for  the  stress  intensity  factor  K^.  It  has  been  found  most  accurate  to 
use  the  tensile  component  of  stress  cr^  with  © = 0 (then  r = x , the 
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distance  along  the  x-axisj.  Apparent  stress  Intensity  factors 
calculated  from  Kq.  with  W = 0, 


are 
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If  the  stresses  from  the  finite  element  analysis  were  exact,  then  a 'onique 
value  of  K-,  would  result,  ’iowever,  extrapolation  is  necessary  since  non- 

sing\ilar  finite  element  analysis  is  not  exact  near  the  crack  tip.  Plotting 

•» 

Kj  as  a function  of  r and  extrapolating  the  linear  portion  of  the 
relationship  to  r = 0 gives  an  estimate  of  the  stress  intensity  factor  at 
the  crack  tip. 

The  extrapolation  can  be  used  only  if  there  is  an  overlapping  range 
of  validity  of  the  singular  solution,  Eq.  (23a),  and  of  the  finite  element 
solution.  An  overlapping  range  of  validity  is  obtained  only  by  extending 
the  range  of  the  numerical  solution  near  the  crack  tip  by  using  a fine  and 
computationally  costly  finite  element  grid. 

E.2  Displacement  matching 

The  displacement  matching  method  is  similar  to  the  stress  matching 
method,  c^orresponding  values  of  displacements  and  position  coordinates 
are  substituted  into  Eq.  (2^4^)  to  solve  for  K^.  It  !ias  been  found  most 
accurate  to  use  the  vertical  component  of  the  displacement  v with 
0 = m so  that 
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Extrapolation  to  r =•  0 gives  an  estimate  of  the  stress  Intensity  factor 
at  the  crack  tip. 
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B.5  Knergy  release 


The  relation  between  linear  elastic  strain  energy  release  rate  G 
and  the  Mode  I stress  intensity  factor  is 


G = 


(29) 


1 T 

where  G = — [f}  [u)  is  the  strain  energy,  A is  the  projected  crack 

area  and 


1.0  , generalized  plane  stress  , 

O = t (30) 

1.0-  V,  plane  strain 

For  "fixed  grips,"  G = - Xi/3a,  with  \s/Mk  < 0 and  the  energy  for 
the  formation  of  new  crack  stirfaces  is  provided  by  the  decrease  in  the 
strain  energy  of  the  body  as  a result  of  crack  propagation.  For  fixed 
boundary  forces,  G = ii/VV,  with  U/  V\  > 0 and  the  energy  for  the 
formation  of  new  crack  surfaces  is  provided  by  the  decrease  in  pot.ential 
energy  of  the  external  loads  due  to  crack  propagation. 

The  energy  release  rate  can  be  evaluated  in  two  ways:  (l)  the 

compliance  difference  [25l ; and  (2  ) the  Bice  .T-integral  [26].  In  the 
compliance  difference  method,  the  release  rate  Ci  is  calculated  by  re- 
placing the  partial  differential  operator,  V in  Eq.  (29),  by  the 
finite  difference  operator  A/.\A  and  measuring  the  strain  energy  of  the 
body  with  two  different  crack  areas.  In  the  Bice  J-integral  method,  the 
release  rate  is  related  to  the  path  independent  J-integral 


UO 


J 


(udy-T*^ds)  = G 


(51) 
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where  u is:  the  strain  energy  density,  T is  the  traction  vector,  u is 
the  displacement  vector,  x and  y are  the  cartesian  coordinates  v’ith 
origin  at  the  crack  tip,  is  the  integral  path  surrounding  the  crack 
tip  and  ds  is  an  element  of  the  arc  along  1'  . 

The  advantage  of  the  energy  release  rate  method  is  that  fine 
finite  element  grids  are  not  required  at  the  crack  tip.  The  Rice  J- 
integral  method  requires  only  one  structural  analysis  to  calculate  G, 
while  the  compliance  method  requires  two  aneilyses. 

B.4  Embedded  singularity 

In  this  method  singular  near  crack  tip  elements  are  used  with  non- 
singular elements  to  model  the  stress  singtilarity  at  the  crack  tip.  Two 
types  of  singular  elements  are  used,  cracked  elements  and  collapsed  ele- 
ments. With  cracked  elements  the  stress  intensity  factor  is  a degree  of 
freedom  of  the  system  and  is  obtained  directly  by  solving  the  equilibritim 
equations.  With  collapsed  elements  the  stress  intensity  factor  is  obtained 
Indirectly  by  displacement  matching. 

For  numericeil  calciilatlons  of  developed  stress  intensity  factors 
in  single  mode  problems,  the  energy  release  rate  methods  are  best.  Stress 
and  displacement  matching  require  fine  finite  element  grids  near  the 
crack  tip,  while  the  embedded  singularity  methods  require  special  sin^rular 
elements  at  the  crack  tip.  For  analysis,  the  release  rate  calculation  by 
the  Kice  J-integral  is  best  since  only  one  stnictural  analysis  is  required. 
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However,  for  design,  the  cmpliance  method  is  best  since  a minimm  of  two 
analyses  are  already  required  for  gradient  calculations,  and  integration 
along  a path  around  the  crack  tip  is  not  required  as  with  the  J-integral. 

C.  Fracture  Constraint  Gradient  Calculations 

Developed  opening  mode  stress  intensity  factors  at  the  crack  tip  are 
calculated  using  the  linear  eleastic  strain  energy  release  rate  method  (the 
compliance  method).  As  was  already  mentioned,  for  munerical  ceilculations 
the  partial  derivative  operator  )/  A in  Eq.  (29)  is  replaced  by  the 

finite  difference  operator  A/\A  and  0 is  obtained  by  taking  the  differ- 
ence in  strain  energy  of  the  body  with  two  different  crack  lengths,  and  = 

+ AA  (the  derivative  applies  to  a mean  crack  length  A = (A^  t-  A,  )/2  ). 
Therefore,  two  complete  structural  analyses  are  needed.  These  two  analy- 
ses also  provide  all  required  constraint  gradient  information. 

Differentiating  Eq.  (29)  with  respect  to  the  design  variables  gives 
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Substltxiting  Eq.  (3^)  into  Eq.  (33)  gives 


^ = r-^  I I 

IdT  40  K,  ' 


pl''’  ' ^i}il 

^ j 


2 A 4M  ) 

4a  ' :a  ' ^dT  / 

I ''  J 


(35; 


At  the  beginning  of  each  redesign  iteration  the  displaceinents , ais- 
placement  gradients  eund  strain  energy  are  calculated  with  crack  area  A - A^. 
'fhe  developed  opening  mode  stress  intensity  factor  is  obtained  by  closing 
the  crack  so  that  A = = A.,  - 'lA,  calculating  the  displacements  and 

strain  energy  and  completing  the  calculation  of  Eq.  i‘P9).  (The  crack  front 
is  closed  by  suppressing  one  or  more  nodes  along  the  crack  line;  the  stiff- 
ness matrix  and  stiffness  gradients  at  A = A^  are  obtained  by  imposing 
the  boundary  conditions  corresponding  to  the  closed  nodes  on  the  stiffness 
matrix  and  stiffness  gradients  at  A = A„ . ) If  the  fracture  constraint  is 

L. 

active,  then  the  displacement  gradients  at  A = A^  are  obtained  and  the 
fracture  constraint  gradient  calcination  of  Eq.  (35)  is  completed.  Grad- 
ients of  active  displacement  and  stress  constraints  are  obtained  only  at 


D.  Numerical  Results 

As  a typical  example,  the  design  of  a center  cracked  stiffened 
aircraft.,  fuselage  panel  is  considered.  The  crack  is  assumed  to  lie  be- 
tween two  frame  members  in  a sjmmetric  manner  so  that  only  a quarter  of 
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the  panel  need  be  considered  (Fig.  f').  A uniform  tensile  stress,  approxi- 
mately equal  to  the  hoop  stress  of  a pressurized  iXt-10  fuselage  panel 
fl7,  1^',  1-]  is  applied. 
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The  finite  element  grid  for  the  skin  jianel  containing  the  crack  in 


sliown  in  Fig.  7 (region  A1?C>  ns  shown  in  Fig.  (i).  A total  of 
finite  elements  is  used  to  model  the  stiffened  panel;  14  two-node  truss 
elements  to  model  the  frsune,  33  truss  elements  to  model  the  longerons  and 
30  planar  three-node  triangles  and  35  planar  Isoparametric  quadrilater- 
als to  model  the  skin  pcinels.  Tiie  15'd  finite  elements  are  linked  to 
three  Independent  design  variables.  The  14  truss  elements  are  linked  to 
(the  frame  area),  the  remaining  33  truss  elements  are  linked  to 
(the  longeron  area)  and  the  105  planar  elements  are  linked  to  (the 

skin  thickness). 

The  frame  and  longeron  material  is  .Aluminiim  2(P4_X3  with  elastic 

/• 

p 

mod’olus  K = 10  X 10  psi*>  Poisson  ratio  v = 0.3,  specific  weight  o = 

•2 

0.1  pci.,  and  yield  stress  = + 37.0  x 10^  psi.  Hie  skin  materieJ.  is 
Aluminum  7073-T^  with  a = + 80.0  x 10^  psi.  The  variation  of  the  crlti- 

y 

cal  opening  mode  stress  intensity  factor  with  thickness  for  Aluminum 

7075-T8  is  shown  in  Fig.  8 [27].  (Secondary  effects  such  as  specimen  width, 
stiffener  spacing,  etc.,  are  neglected.)  For  thicknesses  less  than  approxi- 
mately 0.45  in.  is  constant  and  approximately  equal  to  57.0  x 10^ 
psi. ^ in.  The  plane  strain  critical  opening  mode  stress  Intensity  factor 

3 — ' 

is  taken  to  be  27.0  x 10  psi. v in.  for  thicknesses  greater  than 
approximately  1.0  in. 

The  fail-safe  structural  integrity  requirement  sjjecifies  the  resid- 
ual strength  of  the  cracked  panel  to  be  sufficient  to  arrest  the  crack  and 
permit  damage  detection  in  a "walk  around"  inspection  [20].  The  residual 
strength  or  residual  stress  is  the  applied  stress  required  to  cause  fast 
fracture  of  an  existing  crack  of  a given  length.  The  residual  stress  cr^ 


of  an  unstiffened  panel  decreases  with  increasing  crack  length;  therefore, 


t 

k 


1} 
: ( 


any  crack  which  propagates  rapirily  will  not  be  arrested  and  will  continue 
to  propagate.  The  residual  stress  of  a stiffened  panel  initially  decreases 
with  Increasing  crack  length  but  then  Increases  as  the  crack  approaches 
the  frame  member  since  the  frame  member  relieves  the  stress  at  the  crack 

tip.  Thus,  a fast  propagating  crack  in  a stiffened  panel  may  be  arrested  j 

if  the  frame  members  are  sufficiently  stiff. 

The  increase  in  the  stress  in  the  frame  member  as  the  crack  approaches 
the  frame  may  cause  failure  by  yielding  of  the  frame  material.  Thus,  two 
distinct  failure  mechanisms  are  operative,  fast  fracture  of  the  skin  (the  i 

skin  criterion)  and  yielding  of  the  frame  (the  frame  criterion).  Both  I 

k 

failure  mechanisms  must  be  prevented  in  a fail.- safe  design. 

i 

The  residual  stress  for  the  center  cracked  panel  with  initial  design  | 

o o j 

= O.S  in'),  =0.5  in^)  and  = 0.0?  in.  is  shown  in  Fig.  y.  ;The  j 

initial  design  approximates  test  panels  of  Refs.  17,  l8  and  19.)  The 

residual  stresses  for  both  the  skin  and  the  frame  are  shown,  (■'or  compari- 

son  purposes,  the  residual  stress  for  the  skin  of  sin  unstiffened  panel 

with  skin  thickness  = 0.07  in.  is  also  shown.)  Clearly,  the  frame  failure 

mechanism  will  not  be  active  in  this  initial  design.  i 

) 

(The  residual  stress  curves  were  obtained  by  closing  the  crack  face, 
successively  suppressing  one  node  at  a time  along  the  crack  line,  and  cal-  ' 

culating  the  strain  energy  for  each  crack  length,  developed  stress  inten-  I 

I 

sity  factors  were  compared  to  the  critical  stress  intensity  factor  and  the  ! 

applied  stress  scaled  to  cause  fast  fracture  of  the  crack.  The  finite 
element  grid  is  finer  near  the  frame  than  near  the  center  line  of  the  panel 
because  of  the  increasing  influence  of  the  frame  on  the  approaching  crack. 

U8 
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Fig  9 Resiriudl  Stress  of  Center  Cracked  Panel 
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Thus,  the  residual  stress  calculations  are  more  accurate  near  the  frame 
because  of  the  smaller  crack  length  Increments  used  to  determine  the 
gradients. ) 

llie  residual  stress  usually  attains  a maximtim  value  when  the  crack 
just  proj>agates  past  the  frame  centerline  with  the  (fastened)  frame  re- 
maining intact.  This  situation  cannot  easily  be  simulated  in  the  present 
study  because  truss  elements  are  used  to  model  the  frame  members.  There- 
fore, the  maximum  residual  stress  is  assumed  to  occur  with  crack  half- 
length  a = (lO.O  + 9.75)/2.0  = 9-^75  in.,  the  maxim\mi  crack  length,  that 
can  be  considered.  The  uniformly  applied  tensile  stress  is  scaled  so  ttiat 
the  developed  stress  intensity  factor  with  a = 1.^7b  in.  was  equal  to  the 
critical  value  = 57.0  x 10^  psi.xin.  This  scaling  gives  an  applied 
stress  CT  = 21,70U.4  psi. 

Q. 

A line  of  constant  stress  equal  to  the  applied  stress  intercepts  the 
residual  stress  curve  (Fig.  9)  at  a = 9.^75  in.  and  a = 5.0  in.  Thus, 
a = 5*0  in.  is  the  length  of  the  longest  crack  which  will  be  arrested  at 
the  applied  stress,  and  is  denoted  a^,  the  threshold  crack  length  [2P]. 
Longer  cracks  with  a > a^  will  be  arrested  if  the  applied  stress  is  less 

than  a = 21,704.4  psi. 

a ' 

The  optimum  design  problem  is  stated  as  follows:  minimize  the 

weight  of  the  stiffened  panel  while  maintaining  for  crack  length  a = 

9.875  in.  the  residual  stress  greater  than  or  equal  to  the  applied  stress 
a = 21,704.4  psi.,  or  equivalently,  maintain  with  a = 9* 875  in., 

Kj  > K^.  Gtress  constraints  on  longeron  and  skin  yielding  were  included 
as  well  as  the  stress  constraints  on  frame  yielding.  The  initial  design 
with  =0.5  ini,  =0.5  In^^l, 
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and  Dj  = 0.07  in.  weighed  1(- .9C  lbs. 


Minimum  limits  on  the  lesign  varlaiiles  were  arbitrarily  chosen  (7.  pt-rcent 
of  the  initial  values);  maximum  limits  were  not  imposed. 

Ti\e  final  design  with  0 = 0.688  in'i,  T',,  - II,  - 0.?;i6  in.  »in‘i 

min 

D,  = D-.  = 0.0504  in.  weighed  14.53  lbs.  Tliis  final  design  is  fully 

^ ^min 

converged  with  the  design  factors  (Eq.  9)  ^ 1.0,  < 1.0  and  < 

1.0  and  was  obtained  in  seven  iterations  (requiring  approximately  I6  CPU 
seconds  for  execution  on  the  UCLA  Campus  Computing  Network  IBM  360/9I 
compiling  system  with  single  precision  arithmetic  and  the  FORTFATI  G 
compiler).  Only  the  fracture  constraint  was  active  during  the  design 
process. 

The  frame  relieves  the  stress  at  the  crack  tip  and  reduces  the 
developed  opening  mode  stress  intensity  factor  at  the  crack  tip.  For  the 
initial  design,  “ -O.89I  where  g^  = - l.C  is  the  con- 

straint on  the  developed  stress  intensity  factor.  Increasing  the  skin 
thickness  is  nearly  as  effective  in  decreasing  the  developed  stress  intens- 
ity factor  at  the  crack  tip  as  increasing  the  frame  area;  for  the  initial 
design  = -0.779*  However,  increasing  the  skin  thickness  has  a 

greater  effect  on  weight  than  increasing  the  frame  area  since  W = 
o where  5w/ = 6.40  and  'w/  = 1U8.0  (the  weight 

gradients  are  constant).  For  a given  value  of  the  Lagrange  multiplier  A 
associated  with  the  fracture  constraint,  the  design  factor  =231^  and 
1-  - 4.8  V L,  or  ^L,  - .208  '“'I,.  Thus,  by  Eqs.  10  and  11  , the  skin 

and  the  longerons,  which  have  little  effect  on  the  developed  stress  intens- 
ity factor, are  quickly  reduced  to  minimum  value.  The  minimiim  value  chosen 
for  the  skin  will  have  a large  effect  on  final  weight. 
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The  residioal  stress  cur/e  for  the  final  design  is  shovm  in  Fig.  9. 
Only  the  skin  criterion  is  active.  The  threshold  crack  lerigth  for  the 
final  design  vrith  cr^  = 21,70it.t  psi  is  = ^*.0  in.  Thus,  the  optimi- 
zation decreases  the  panel  weight  (approximately  ly.5  percent)  but 
also  decreases  the  threshold  crack  length,  the  length  of  crack  which  will 
just  be  arrested  at  the  applied  load.  Durability  requirements  of  the  type 
proposed  in  Ref.  20  can  also  be  included  in  the  present  opt.imization 
algorithm.  For  details  see  Appendix  D. 

E.  Disc\ission 


The  results  for  the  center  cracked  stiffened  panel  show  that  sig- 
nificant weight  decreases  are  possible  while  maintaining  the  structural 
integrity  of  the  panel.  Similar  results  were  obtained  for  a center  hole 
cracked  stiffened  panel  (modeled  by  deleting  the  longeron  four  inches 


above  the  crack  line  of  the  center  crack  model  and  inserting  an  eight  inch 
diameter  hole).  Again,  the  weight  decrease  is  obtained  at  the  cost  of  an  f; 


increased  sensitivity  to  initial  crack  sizes  (i.e.,  a decrease  in  thresh- 


old crack  lengths 


a 

a 


The  criteria  for  fracture  resistant  aircraft  structures  proposed 


in  Ref.  20  includes  a dxirability  requirement  as  well  as  a structural  inte- 
grity requirement.  The  structural  integrity  requi."ement  ensures  that 
cracks  will  not  grow  to  critical  size  during  the  life  span  of  the  stricture. 
For  example,  the  durability  requirement  proposed  in  Ref.  20  for  fail-safe 
design  states  that  an  initial  crack  of  half  length  a = O.T’R  in.  not  grow 
to  critical  size  in  two  life  spans  of  the  aircraft.  The  growth  rate  is 
calculated  ignoring  structural  design  details.  This  requirement  ensures 
against  '’xpenslve  repair  of  a panel  failure  durlrig  the  life  span  of  t>;e 
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Iho  aircraft..  i)ural  ilLty  requirement  constraints  can  i>e  included  in  tfie 
present  metho<l  at  some  aflditionixl  cotnputal  ional  cost. 

Hie  results  oP  the  present  study  were  oV.tained  by  modeling  the 
rrane  and  longerons  only  by  truss  members.  More  accurate  results  coiild 
be  obtained  by  a more  detailed  mcnieling  of  the  frame  and  longerons  and 
their  cormections  to  the  sheet.  rmproved  modelirif.-  v’o\ild  permit  study  of 
crack  stoppers  and  frame  miCmbers  of  different  materials  and  a].lcn»-  consid- 
eration o'"  secondary  effects  such  as  st>~ir4;er  eccentricity,  river  spaciiig, 
etc.  h , --.’fi}.  "This  more  detailed  (anci  more  costly/'  moJeliryi;  would  allow 
consideration  of  more  practical  cases  wf’.ere  cracks  span  two  skin  panels, 
since  cracks  usually  do  not  start  at  the  middle  of  a panel  hut,  rather, 
start  at  stiffener  attachments  and  mm  into  one  ov  both  adjacent  skin 
panels . 

The  present  results  were  obtained  with  the  total  mimber  of  finite 
elements  linked  to  three  independent  design  varLal  les.  Improved  results 
may  be  expectedi  if  a less  severe  linking  were  imposed  and  a larger  n'jmber 
of  design  variables  were  considered,  especially  near  proliable  crack 
locations . 

The  designs  obtained  by  the  present  method  are  > ased  on  the  exist- 
ence of  a crack  at  a specified  location  and  with  a si/ecifi^d  geometry  and 
assumed  opening  tno<le  tip  deformation.  Several  cracks  at  different  loca- 
tions should  be  post\ilated  for  a more  realistic  application  of  this  method. 
A finite  element  analysis  for  each  crack  could  he  implemented  and  included 
in  the  present  design  algorithm.  Tlius,  the  metho<i  can  be  applied  to  a 
wide  range  of  design  problems  reqiiiring  the  inclusion  of  fract'ure 
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constraints. 

Finally,  it  should  he  noted  that  for  this  problem  the  fractixre 
mechanics  constraint  was  the  dominant  design  consideration,  i.e.,  stress  auid 
displacement  constraints  were  not  active.  For  other  applications  it  is 
clear  that  stress,  displacement  and  side  constraints  may  become  active, 
but  clearly,  design  for  fracture  prevention  is  a very  important  effect  and 
must  be  incorporated  in  automated  structural  design. 


V.  CLOyi/RK 


The  primary  research  objective  In  the  study  of  optimality  criteria 
approach  to  automated  structural  deEi^n  has  been  to  develop  an  efficient 
design  algorithm  capable  of  treating  large  structural  systems  witr.  several 
hundred  or  more  design  variables.  In  addition,  the  algoritlim  has  been 
developed  with  a capability  for  dealing;  with  mlnim’jm  weight  design  prot - 
lems  involving  stress,  displacement,  side  constraints  on  otner  desiy,n 
variables  and  constraints  against  fracture. 

As  described  in  this  Interim  Report  these  objectives  have  been 
attained.  The  present  algorithm  is  capable  of  determining  very  efficient 
designs  of  a wide  range  of  aerospace  structurfts,  subjected  to  miltiple 
static  l(5ad  systems  and  constrained  against  several  forms  of  inadmissatle 
behavior.  The  algorithm  has  been  tested  with  severs''  well  known  difficult 
elastic-truss  design  problems  and  has  al.so  been  exercised  on  design  prob- 
lems involving  built  up  assemblages  of  two-dimensional  stiffened  panel 
structures. 

The  present  algorithm  has  been  shown  to  produce  near  optimum 
designs  with  a significantly  smaller  computer  effort  than  the  best  avail- 
able mathematical  programming  techniques.  Although  the  present  alg.orithra 
does  require  more  computational  effort  than  other  heuristic  "optimality” 
criteria  methods  (since  the  algorithm  is  based  on  t’ne  trrie  mathematical 
optimality  criteria,  the  KulJi-Tucker  test),  the  designs  produced  by  the 
algorithm  are  at  a local  optira'jm  point.  This  capability  has  been  shown 
fll  to  produce  designs  v;hich  are  s\q)erior  to  designs  developed  using  the 
heuristic  optimality  criteria  methods. 


Perhaps  the  most  important  aspect  of  the  research  reported  herein 
is  the  development  of  an  automated  structural  design  capability  for  treat- 
ing constraints  against  fast  fraot\ire  due  to  praj>agation  of  existing  struc- 
tural cracks.  As  indicated  by  the  results  in  Section  IV,  the  inclusion  of 
fracture  mechanics  constraints  can  have  a dominant  effect  on  optimum 
structural  design  since  the  structures  designed  with  this  constraint  in- 
cluded were  not  constrained  by  the  more  conventional  naximtim  and/or  minimum 
limits  on  allowable  displacements  and  stresses,  ^'ased  on  these  results,  it 
is  recommended  that  research  continue  in  the  area  of  automated  design  of 
aircraft  structures  for  fracture  resistance  and  safe-life  requi'-ements . 
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APPENDI-X  A 
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A computer  program  for  efficient  minimum  weigf.t  design  of  large 
scale  structural  systems  is  described  in  this  section.  Tlie  prograr* 
utilizes  finite  elements  for  stress  analysis  and  the  huhn- Tucker 
necessary  conditions  for  design. 

Five  types  of  finite  elements  are  applied  to  the  stress  analysis. 
These  are : 

(l ) the  three-dimensional,  two-node  truss  element  (TP.USSf  ); 

{'d  ) the  three-dimensional,  three-notle  subparametric  truss 
element  (THUSS:^); 

(:5)  the  planar  isotropic  three-node  constant  strain  triangle 
(CFT); 

{‘0  the  planar  isotropic  four-node  isoparametric  quaurilateral 
(g,UAl)4);  and 

(5)  the  planar  isotropic  eight-node  isoparametric  quadrilateral 
(QUAD8). 

!X*tails  for  the  stiffness  matrix  calculations  are  given  in  A^jpendix 

Four  types  of  behavioral  constraints  are  included  in  the  present 
version  of  the  design  program.  These  are: 

(1)  limits  on  the  tensile  and  compressive  stresses  in  the  truss 
elements  and  limits  on  the  von  flises  equivalent  stresses  in 
the  planar  elements ; 

('<; ) limits  on  nodal  displace  nents ; 

(d)  limits  on  relative  nodal  displacements;  and 

(><■)  limits  on  developed  opening  mode  stress  intensity  factors. 
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In  addition,  !riinimum  nna  mnxim\im  side  limits  on  the  design  varia;  les  are 
included . 

Tlie  program  ha.s  91  subroutines  including  the  main  subroutine.  '11 
routines  are  written  in  FORTRAIi  IV.  Only  internal  storage  is  usea.  Ail 
arrays  are  included  in  a single  blank  common;  the  length  of  each  a''ra2'  is 
dynamically  allocated  in  terms  of  the  pertinent  input  paramet'^rs . In 
addition,  there  are  b2  labeled  common  blocks.  All  real  variables  and 
constants  are  declared  as  single  precision  numbers. 

Program  Structure 

The  structure  of  the  program  is  illustrated  in  Fig.  A1 . .'ubrout.ine 
MAIN  is  listed  in  Fig.  A2 . Overall  control  of  the  arislysis  and  desigr.  is 
accomplished  in  fWIK.  AJ.1  control  transfers  are  retiirned  to  .’AMT'i. 

The  four  ma.ior  functions  of  the  program  are;  (l)  (iata  input;  ( ) 
analysis;  (5)  constraint  evaluation;  and  (-i)  redesign.  The  first  major 
function,  data  input,  is  accomplished  in  the  steering  subroutine  lATAIIi. 
Besides  data  input,  bATAIU  accomplishes  the  dynamic  allocation  of  all 
arrays  and  the  calculation  of  the  stiffness  gradients  (and  storage  in 
compacted  form ) . 

The  second  major  function  is  a finite  element  stress  analysis, 
including,  when  sp>ecified,  the  calculation  of  developed  opening  mode  stress 
intensity  factors  by  the  compliance  method.  The  system  stiffness  matrix 
is  synthesized  in  STIFF,  the  system  displacements  calculated  in  GLIl^'F 
and  the  member  stresses  calculated  in  STRESS.  If  <'"acturo  constraints  are 
active  ( IACTIV(U  )=1 ) then  two  analyses  R’-e  required,  one  analysis  with 
crack  area  A.^  and  the  other  analysis  with  crack  area  = ■V  - A.A. 
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Bounciarj'-  conditions  correspond Ing  to  suppressed  nodes  (to  close  the  crack 
from  A.,  to  ) are  imposed  in  (’BOUND.  Otress  intensity  factors  are 
then  calculated  in  CAK’KI. 

The  third  major  function  is  the  evaluation  of  the  constraints, 
including  the  identification  of  near  active  constrains  in  ORDER.  Gradients 
of  the  near  active  constraints  are  then  obtained  in  GPDEIJT  and  dependent 
constraints  are  eliminated  in  ELEiilT. 

The  fourth  major  function  is  the  iterative  redesign.  The  calcula- 
tion of  the  Lagrange  m\iltlpliers  and  the  design  factors,  the  imposition  of 
move  limits,  the  trial  resizing  and  the  elimination  of  Inactive  constraints 
are  accomplished  in  steering  routine  liEDESN. 

input  Data  Preparation 

Required  input  data  cards  are  described  fully  in  Appendix  B.  Tlie 
following  description  is  only  for  that  portion  of  data  which  differs  from 
the  data  required  in  ordinary  finite  element  programs. 

i . Design  variable  linking 

Element  sizing  variables  A^  (the  areas  of  the  truss  elements  and 
the  thicknesses  of  the  planar  elements)  are  linked  to  design  variables 
through  the  linking  matrix  a^ . so  that 


A 


i 


ii 


1,2,...,.. 
1,2, ...,m 


(Al) 


?)enernlized  linking  is  provided i that  is,  different  element  types  may  be 
linked  to  the  sane  design  variable  and  elements  may  be  linked  to  one  or 


ti6 


morp  fiesign  variables.  For  eacli  element  it  is  necessary  to  specify  the 
cor respond ing  design  varlalile  and  the  value  of  the  linking  coefficient. 


' • Constraint  criticality  buffer 

A constraint  is  strictly  critical  if  g^^  = 0.0  or  if  * 
1.0  = 1.0.  For  numerical  piu-poses,  a constraint  is  assumed  active  and 
governing  the  redesign  if 


FILTER  = 1.0  - P < < 1.0 


(AO) 


where  P is  the  constraint  criticality  buffer  and  FILTER  is  the  variable 
name  of  the  quantity  1.0 -p.  !)efault  values  of  FILTER  (and  the  constraint 
criticality  buffer  p)  are  obtained  by  reading  in  a zero  initial  value. 

5*  Impendent  constraint  filter 

The  linear  equations  to  solve  for  the  Lagrange  multipliers  of  active 
constraints,  Eq.  (l8),  are  functions  of  the  constraint  gradients.  The 
constraints  for  problems  with  symmetry  linking  are  dependent,  and  therefore, 
so  are  the  constraint  gradients  and  the  equations  for  the  Lagrance  multi- 
pliers. Dependent  equations  are  eliminated  by  comparing  constraint  grad- 
ients term  by  term.  Constraint  gradients  are  assumed  equal  if  the  differ- 
ence is  less  than  ESIF  where 


ESIF 


10 


-NS  IF 


(A5) 


and  NSIF  is  the  expected  accuracy  of  the  computer  in  terms  of  the  number 
of  significant  digits. 


b 

I 


j 


I 

V 
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APPENDIX  B 


TNRJT  DATA  DESCHTPTION 


l)  Control  dn'ta,  two  cards. 

Card  1:  (P0A4);  heading 

CardC: 

colurans  1-5;  number  of  nodes, 
columns  6-10;  number  of  element  types. 
colvDiins  11-15;  number  of  elements, 
columns  l6-C0;  number  of  design  variables, 
columns  Pl-25;  number  of  load  conditions, 
columns  26-50;  number  of  design  cycles. 

C)  NcKie  point  data  inpvit,  one  card  per  node  15E15.6,  515). 


columns  1-15 
columns  ll-26 
columns  27-59 
coliimns  40-14 
columns  45-49 
columns  50-54 


coorainate 
Xr,  coordinate. 


coordinate. 


code  for 


code  for  x,. 


code  for  x. 


0,  free. 


1,  constrained. 


5 ) Element  input . 

one  card  per  element  type  (4ib). 


columns  1-  ^ 


columns  6-10 


columns  U-IB 


columns  16-20 


element  cards. 


type  number  of  element, 
number  of  elements  of  type. 


blank. 


blank. 


Card  1:  (5E15.6,  15). 


columns  1-15;  elastic  modulus. 
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U) 


columns 

ld-2f;; 

Poisson 

columns 

27-39; 

we ight 

columns 

40-44 ; 

aneLLys  i 

one  card  per 

element 

(815) 

columns 

1-  5; 

node  I. 

colimfins 

6-10; 

node  J. 

columns 

U-15 ; 

node  K. 

columns 

16-20; 

node  L. 

columns 

21-25; 

node  M. 

columns 

26-30; 

node  N. 

columns 

31-35; 

node  0. 

columns 

36-40; 

node  P. 

nodes  numbered  ccrunterclocKwlse 
from  I - P 


two  cards  per  element 

Card  1:  (l^);  design  variable  number. 

Card  2:  (E15.6);  linking  coefficient 
Design  variable  input. 

Card  1:  (e15.6);  default  design  variable. 
Card  P : (EE16.6). 

columns  1-13;  default  minimum  val\ie. 
oolximns  lU-'d6;  default  maximum  value, 
one  card  per  design  variable  (3E13.f^). 


columns  1-13 
columns  1^-Pf 
coliimns  i’7-39 
Load  data  Injrut . 

on'*  canl  per  loaded  nofie  3E16.6) 


design  variable  value, 
minimum  value, 
maximum  valiie. 


Jk 


f) 


coluinns  1-  5j  node, 
columns  6-10;  load  condition, 
columns  11-23;  component  of  load, 

columns  2U-26;  Xo  component  of  load, 
columns  37-^9;  x-  component  of  load. 

Constraint  data  input. 

control  card:  (415);  constraint  activity.] 

columns  1-  5;  stress  activity 
columns  6-10;  displacement  activity 

- relative  displacement 

columns  11-15 ; 


activity 

columns  l6-2C;  fracture  activity 
element  cards  (if  lACTIV(l)=l) . 

Card  1:  (’E13.6). 

columns  1-13;  default  tensile  yield, 
colximns  l4-?6;  defaiili.  compressive  yield, 
one  card  per  element  (2E13.6). 

colimins  1-13;  tensile  yield  stress, 
columns  14-26;  compressive  yielu  stress, 
node  cards  (if  IA0TIV(2 )=l). 

Card  1:  (2K13.6). 

columns  1-13;  default  upper  displacement  limit 
columjis  14-26;  default  lower  displacement  limit 
one  card  per  node  (2E13.6). 

columns  1-13;  upj>er  displacement  limit, 
columns  l4-06;  lower  displacement  limit. 


0,  nonactive, 

1,  active. 
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relative  displacement  cards  (if  L\CTTV(.5)=l). 

Card  1:  relative  displacement  limit 

one  card  per  constraint  (CiIO). 
columns  1-  5>  node  A. 

columns  6-10;  node  B. 

columns  11-15;  load  condition, 
columns  16-00;  direction  code  j 1 

columns  01-25;  Xr,  direction  code 

columns  26-30;  x..  direction  code  | 1 

fracture  cards  (if  IACTIV(4)-=l). 

Card  1:  (E13.6);  intensity  factor  limit. 

one  card  per  load  condition:  (I5);  load  condition  activity. 

7 ) Crack  data  input  ( 1 f IACTIV(4  )=1 ) . 

Card  1 : ( 3K15 • 6 ) . 

columns  1-13;  crack  area  factor, 
columns  14-26;  crack  length, 
columns  27-39;  crack  ext.ension. 
one  card  per  constrained  crack  node  (5I5)« 
columns  1-  5>  botinding  element  tjT)e. 

columns  6-10;  bounding  element  number, 

columns  11-15;  design  variable  number, 
columns  16-20;  initial  crack  front  node, 
columns  21-25;  final  crack  front  node. 

6)  Constraint  elimination  data  input  (5I5»3K13»6). 

coliirans  1-5;  number  of  constraint  gradients  calculated. 


relative  displacement 
not  active. 


relative  displacement 
active . 
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columns  blank;. 


columns  11- lb;  elimination  option, 
columns  lf^-?:R;  Kf-D’. 
columns  blank, 

columns  FTLTKR. 

9)  Move  limit  data  i?Vi ,2E1'^,2,) . 

columns  1-  b ; blank. 

columns  ^<-10;  fully  stressing  option. 

columns  11 -R 5;  move  limit  A. 

columns  5^-;  fiiUy  stressing  move  limit  I,. 

10)  Print  switch  data  input  (5T5) 

columns  1-  5i  member  area  print  switch, 
columns  R-10;  nodal  displacement  print  switch, 
columns  11-lb;  member  stress  prljit  switch, 
columns  16-20;  constraint  print  switch, 
columns  21-2b;  gradient  print  switch. 

switch  = 0,  no  print, 
switch  = 1,  print. 
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A pram)  IX  c 


FINITE  ELEMENTS  USED  DJ  STRESS  AILUYSIS 

The  stiffness  relations  for  a linear  elastic  structure  are 

fK]fu]  - fF]  (Cl) 

where  the  system  stifPness  matrix  [K]  is  obtained  by  direct  addition  of 
the  element  stiffness  matrices  with 

[K]  = I tn]^  lE][B]dv  . (ra; 

- V 

In  F,q.  (C2  ),  [b]  is  the  linear  strain  displacement  matrix  for  element  i, 

I El  is  the  matrix  of  elastic  constants  and  V is  the  volume. 

Brief  derivations  of  the  stiffness  properties  of  the  finite  elements 
used  in  the  present  version  of  the  design  program  are  given  belcw.  The 
two-node  truss  element  (TRUES? ) is  a special  case  of  the  TRUSS?  element 
and  the  three-node  triangle  (CST)  and  the  four-node  quadrilateral 
(QUAI^)  are  special  cases  of  the  Q,UAr»8  elements;  therefore,  the  TRUSS?, 
the  CST  and  the  QUAUL  elements  are  not  described  separately.  Idirther 
details  of  the  stiffness  calculation  can  be  found  in  many  texts,  including 
Ref.  ?9. 

Three-Node  Three-Dimensional  Truss  (TRUSS?) . 

The  three-node  subparametric  truss  element  is  shown  in  Fig.  Cl.  The 
coordinates  x^,  x,  and  x^  are  system  cartesian  coordinates,  S is  a 
local  coordinate  along  the  truss  varying;  from  0 to  L (the  truss  length) 
and  r is  a natural  corrdinate  along  the  truss  varying  from  -1  to  »1. 


7? 


T>ie  d ispIacRment 
natural  corrdinnte  r, 


u alonf^  the  truss  is  assumed  qiiadrai.ic  in  the 
that  is, 


u = 


T (l  - r ) - (l  - /'  ) |u. 


|^(l+rJ-^(l-r‘')'lu^j  (l 


r^  )u . 


(C3J 


where  u^,  u^  and  u^  are  the  system  displacements  at  nodes  i,  j and 


k,  respectively.  The  system  nodal  displacements  may  he  vm-itten  in  terms 

of  the  components  along  the  system  cartesian  coordinates 

etc.)  and  the  projected  lengths  of  truss  along  the  system  coordinates 


(L^,  I+,  and 

where 

L = 

O 

x/  ^ 

+ 1?  + 

);  that  is, 

^i 

"l 

u.., 

li- 

h 

""il 

L " 

L 

L 

= "J1 

^.12 

u. 

r. 

T 

"J5 

L 

(Cu) 

= \l 

b. 

L • 

The  linear  strain  displacement 

relation  is 

£ 

du 

dS  ' 

\C^) 

The  displacement  is  given  in  terms  of  the  natural  coordinate  r and  the 
strain  displacement  relation  is  in  terms  of  the  local  coordinate  G. 
Therefore,  it  is  necessary  to  relate  tne  natural  coordinate  derivative 
and  the  local  coordinate  derivative  to  evaluate  the  strain.  Py  the 


. 1 


1 


1 

1 


n 

t1 
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chain  rule 


d _ ^ s T — 

dr  dr  dS  dS 


ic6) 


where  J is  the  Jacobian  operator  relating  the  natural  and  local 
derivatives. 

The  local  position  coordinate  S is  linear  in  terms  of  the 
natural  coordinate  r,  that  is, 


S = I a + r)  (C7) 

and  J = dS/dr  = L/''’ . The  order  of  the  position  coordinate  in  r is 
less  than  the  order  of  the  displacement  function;  hence,  the  term 
subparametric 

The  strain  is 


dS 

— t" 

■1  ^ 

— 0 

dr 

o 

e. 

L 

r + ( J + 

r )u^  - 2r  uj  . 

(C8) 

Substituting  for 

Uf, 

u and 

from  Eq . 

('y.  } gives 

e 

(cy) 

where 


For  constant  truss  area  A, 


r 


IB]^  E[B]ds 
[B]”^  E[B]  k d.r  . 

I . 


(C15) 


ExiJlicit  integration  is  replaced  by  two-point  Guass  quadrature  so 

that 


[K] 


i 


A E L 
2 


I 


W^[B]^  [P] 

J li 


(Clt) 


where  W , is  the  weightiTig  function  at  integration  point  j and  [b1  . 

J J 

is  the  strain  displacement  matrix  evaluated  at  integration  point  j. 
Eight-Node  FLanar  Quadrilateral  (QUAD8). 


The  eight -node  isoparametric  planar  element  is  shown  in  Fig.  C’ . 
The  coordinates  and  x,  are  system  cartesian  coordinates  and  r 
and  s are  natural  coordinates  varying  from  -1  to  +1, 


The  geometry  of  the  element  in  the  cartesian  coordinates  is 
"mapped"  into  a simple  square  in  natural  coordinates  for  ease  in  inte- 
gration for  the  stiffness  coefficients.  The  mapping  or  interpolating 
functions  relate  the  cartesian  and  natural  coordinate  positions  so  that, 


’'l  = 

fit,} 


(CIB) 
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where 


I 


^ (1  - r)(l  - s)  - i (1  - r)(l  - s'  ) - ^ (1  - r"  )(l  - s) 

^ (l+r)(l-s)  - ^ (l^r)(l-s‘') 

I (If  r)(l+  s)  - J (1-  /)(lf  e)  - ^ (if  r)a-  s^) 

r (l-r)(lfs)  - f (l-r^)(lfs)  - i (l-r}(l-s"') 

4 4 4 

I (l-r^)d-s) 

^ (l f r)(l  - ) 

I (1  - r"  )(lf  s) 

I (1  - r)(l  - ) 


(C16) 


and  fx]  Is  the  vector  of  nodal  x^ 

The  displacements  and  Up, 

terms  of  the  natural  coordinates  are 


coordinate  positions. 

along  the  x^  and  x^  axis  in 


'^l 

. 

[0]'"  1 

a. 

(of 

fa’f  _ 

fd,  I 


(C17) 


where  fctl  are  the  interpolating  or  mapping  functions  and  fu^')  and 

I tu  1 are  the  vectors  of  the  system  nodal  displacements  along  the  x^ 

ami  X,  axes,  respectively.  Tlie  geometry  and  the  deformation  are  the 

same  order  in  the  natural  coordinates,  that  is  fc,}  = fc;  Hius,  the 

g 
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(-'Ji 


(ci8) 


or 


= 


E 

11 

(0,^- 

(oi’‘ 

1 

1 

- 

1 

9-^ 

fu,' 

* c.  ' 

= [B]fu] 


(C19) 


It  is  necessary  to  relate  the  natural  coordinate  derivatives  and 
the  system  coordinate  derivatives  to  evaluate  the  strains.  By  the 
chain  rule 


^ I 

1 

iXi  V 1 

^ S 1 

V 1 

= [J] 


(Chi) 
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1 


1 


Therefore, 


1 ~)  » 

1 1 

i; 

= (J]‘^ 

1 i 1 

^ \ 

1 ^ I 

I 

Equation  (C19)  can  be  used  to  obtain  the  strain  displacement  matrix  [B] 
in  terras  of  the  natural  coordinates. 

The  stiffness  matrix  is 


:k]. 


[B]  [E][B]dV  . 


(C^2  ) 


For  constant  thickness  t, 


[Kl.  = t 


[b]  lE][B]dA 


A 

^ r 
1 1 


= t ^ [ [B]”^  [E]  1b1  I J|dr  ds  . 


(F-'i) 


The  explicit  integration  in  the  natural  coordinates  is  replaced  by 
lauss  quadrature  so  that 


I Kl  . = t 


A n 

L L ^ fb]  [B(r^  , )1 

,1=1  1=1 


Be 
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EFFKCT  OF  DESIGN  VARIABLES  ON  LIFE  OF  STRUCTURES  SUBJECTED 
TO  CONSTANT  AMPLITO])E  CYCLIC  LOADINCiS 

When  linear  elastic  fracture  mechanics  is  applicable  the  stable 

crack  growth  i-ate  cia/dN  during  constant  amplitude  cyclic  loading 

(between  o and  a . ) can  be  expressed  as 
max  min  ‘ 


^ > K ) 

dN  mean 


where  K is  the  number  of  cycles  of  loading,  IK  = Kj  ” “I  ’ 

max  min 

K = ^ (K-r  * ).  The  parameter  K^.  is  the  stress  intensity 

max  min  max 

factor  at  a = cr  and  is  the  stress  intensity  factor  at  cr  = 

max  I , 

min 

o (K^  =0  if  O'  , < O).  Experimental  results  exhibit  a linear 

min  I„.„  min 

relation  between  log  (da/dN)  and  log  (Ak)  for  constant  values  of  mean 
stress  (with  significant  deviation  only  near  final  failure)  [50],  that  is, 


i = 


where  C^  and  n are  empirical  constants  determined  from  the  experimental 
results  with  n usually  in  the  range  2.i  < n < 4.0. 

Solving  Eq.  (fC)  for  dN  and  integrating  between  the  initial  crack 
size  a^  and  the  critical  crack  size  a^  gives  the  number  of  cycles  to 
failure  (the  life  of  the  structure), 

N - ^ . (D3) 

^ '^0  '-  a.  (Ak)”^ 

The  gradient  of  the  life  with  respect  to  a design  variable  is 
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r 


c 


a 

f " 


c 


1 , 


(5)-<  ) 


The  durability  requirement  in  Ref.  20  states  that  a^  = O.l/!^  in. 


P'or  fail-safe  design  of  stiffened  panels  a^  = a^,  the  threshold  crack 


size.  The  stress  intensity  factor  Ak  and  the  gradient  ^K/^D.  can 


be  calculated  as  discussed  in  Section  IV.  Then  the  integral  equation 
for  the  life,  Eq.  (D5),  and  the  integral  expression  for  the  life  gradient, 
Eq.  (dU),  can  be  evaluated  numerically.  To  complete  the  life  gradient 


calculation  it  is  necessary  to  evaluate  da  /dD.. 

c'  j 


In  the  initial  stages  of  crack  growth,  the  crack  length  is  small  and 
the  stiffener  has  little  effect  on  the  near  crack  tip  stresses.  Tlaerefore, 
an  approximate  calculation  for  the  critical  crack  gradient  can  be  obtained 
using  the  analytic  solution  for  a crack  in  an  infinite  sheet. 


cr^ 


TTa 


(n5; 


Equation  (DS)  is  used  only  to  estimate  the  gradient  da  /dD.  in  Eq.  (lk+). 

^ 0 


At  the  onset  of  fast  fracture  a = and  Substituting 


a and  into  Eq.  (D''W  and  differentiating  with  respect  to 

c i C 


Dj  gives 


dK„ 


d 


■'/  Ta  I + 
c 


W I 


TTa 


dj 


Solving  for  da  / dD . giver 
^ J 


1 


* 
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(17) 


?a  2 n/tte  ' 'K,, 

c c C 

a no 


For  the  example  problem  with  design  variables  fp  and 

(and  assuming  the  crack  in  the  stiffened  panel  acts  as  a crack  in  an 


infinite  unstiffened  sheet), 


c = ^ 
^3 


where  q is  the  far  field  line  load.  Therefore, 


0 ^ J 3 

3o 

^ -a/h.,  j = 3 


so  that 


J3  D, 


(Dio) 


where  is  the  Kronecker  delta.  This  result  implies  that  is 

constant  over  a sufficient  range  of  (see  Fig.  8). 

Substituting  Eq.  (DIO)  into  Eq.  (D^),  gives 


"o  I AK^(a^)  ^3 


3^  , 
^ da 

J 


(nil) 


Numerical  results  for  the  stiffened  panel  with  = 0.^  in.,  II,  = 

O 

0.5  in'i  and  = 0,07  in.  and  a = a^  = 5*0  in.  gives 
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showing,  as  expectol,  that  a decrease  in  the  design  variables  will  decrease 
the  life  of  the  structure. 


I 

END  9-77 
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